Utilizaremos los datos “engel” del paquete “quantreg”, datos utilizados en Koenker y Bassett (1982). Se trata de un conjunto de datos que consta de 235 observaciones sobre la renta y el gasto en alimentos de los hogares de clase trabajadora belgas. Se desea obtener un modelo que explique el gasto en alimentos en función de la renta familiar.
En lugar de ajustar todos los modelos robustos que hemos visto, recomiendo utilizar la función “ltsReg” del paquete “rrcov” que permite obtener un resumen del modelo comparable al de la función “lm”.
Comparar los resultados de los modelos lm y ltsReg.
library(quantreg)
data("engel")
head(engel)
## income foodexp
## 1 420.1577 255.8394
## 2 541.4117 310.9587
## 3 901.1575 485.6800
## 4 639.0802 402.9974
## 5 750.8756 495.5608
## 6 945.7989 633.7978
attach(engel)
boxplot(engel)
# Se aprecia que no siguen normalidad, y existen outliers en ambas variables.
shapiro.test(engel$income)
##
## Shapiro-Wilk normality test
##
## data: engel$income
## W = 0.79204, p-value < 2.2e-16
# p < 0.05. Se rechaza la hipótesis nula de existencia de normalidad.
shapiro.test(engel$foodexp)
##
## Shapiro-Wilk normality test
##
## data: engel$foodexp
## W = 0.87471, p-value = 5.45e-13
# p < 0.05. Se rechaza la hipótesis nula de existencia de normalidad.
plot(income, foodexp, pch=20)
abline(coef(lm(foodexp~income)), col='green')
En principio sí parece haber linealidad.
library(WRS)
pbcor(income, foodexp)
## $cor
## [1] 0.9342118
##
## $test
## [1] 39.9758
##
## $siglevel
## [1] 0
##
## $n
## [1] 235
wincor(income, foodexp)
## $cor
## [1] 0.9273018
##
## $cov
## [1] 34445.98
##
## $siglevel
## [1] 0
##
## $n
## [1] 235
La correlación entre ambas variables es muy fuerte y positiva. p-valor menor a 0.05 por lo que rechazamos la hipótesis nula (correlación es 0). Existe una relación significativa entre ambas variables.
# Alimentos (VD) en función de la renta (VI). Es decir, foodexp en función de income.
summary(fitLS <- lm(foodexp~income))
##
## Call:
## lm(formula = foodexp ~ income)
##
## Residuals:
## Min 1Q Median 3Q Max
## -725.70 -60.24 -4.32 53.41 515.77
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 147.47539 15.95708 9.242 <2e-16 ***
## income 0.48518 0.01437 33.772 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 114.1 on 233 degrees of freedom
## Multiple R-squared: 0.8304, Adjusted R-squared: 0.8296
## F-statistic: 1141 on 1 and 233 DF, p-value: < 2.2e-16
Como el intercepto (beta 0) e income(beta1), poseen un p-valor < 0.05 rechazamos la hipótesis nula por lo que ambas son significativos, deben aparecer en la ecuación.
La ecuación quedaría como: foodexp = 0.48518 * income + 147.47539
F-statistic (ANOVA), el p-value es < 0.05 nos aporta que el modelo es significativo.
Los parámetros del modelo intercepto y pendiente son ambos significativos (p-valor < 0.001)
Se nos muestra en F-statistic del paso anterior. Seria: F(1,233)=1141, p-value < 0.001
La bondad de ajuste viene dada por R2 que en este caso es de 0.8304, el modelo es muy eficaz (entre 0 y 1).
# graficamos...
plot(foodexp~income, pch=20)
abline(fitLS)
La linealidad se puede apreciar en el gráfico de los datos anterior. Existe linealidad.
# Gráfico genérico del modelo:
par(mfrow=c(2,2))
plot(fitLS)
library(lmtest)
dwtest(fitLS)
##
## Durbin-Watson test
##
## data: fitLS
## DW = 1.4108, p-value = 2.677e-06
## alternative hypothesis: true autocorrelation is greater than 0
Se usó prueba de Durbin-Test. Se obtuvo: Autocorrelación positiva dwtest < 2.
# Residuals vs. Fitted.
par(mfrow=c(1,1))
plot(fitLS, which=1)
# No aparece ningún tipo de patrón. Parece media 0. No hay heterocesdasticidad ni problemas de linealidad. No hay correlacionalidad en los residuos.
# Se aprecian 3 outliers 59, 105 y 138.
# Normal Q-Q
plot(fitLS, which=2)
# A excepción de los outliers comentados anteriormente, parece que se ajustan los datos a la recta.
# La hipótesis de que los errores siguen una distribución normal parece cumplirse.
# Scale - Location
plot(fitLS, which=3)
# Como en el caso del gráfico primero, no aparece ningún tipo de patrón.
# Residuals vs. Leverage
plot(fitLS, which=4)
# Muestra los 3 outliers que se deben estudiar.
# graficamos distancias de cook
cooks <- cooks.distance(fitLS)
hat <- lm.influence(fitLS)$hat
par(mfrow = c(2,1))
plot(income,foodexp, col=ifelse(cooks > quantile(cooks, .90), 'red', 'black'), pch=20)
plot(income,foodexp, col=ifelse(hat > quantile(hat, .90), 'green', 'black'), pch=20)
par(mfrow = c(2,1))
plot(income,predict(fitLS), col=ifelse(cooks > quantile(cooks, .90), 'red', 'black'), pch=20)
plot(income,predict(fitLS), col=ifelse(hat > quantile(hat, .90), 'green', 'black'), pch=20)
# Parece que hay outliers que influyen en la recta. Se debería trabajar mediante el uso de técnicas robustas (Ya comentado al inicio).
# Uso de validación cruzada para ver el ajuste del modelo.
library(caret)
datos <- engel
train(foodexp~income, data = datos, method = "lm")
## Linear Regression
##
## 235 samples
## 1 predictor
##
## No pre-processing
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 235, 235, 235, 235, 235, 235, ...
## Resampling results:
##
## RMSE Rsquared
## 124.3693 0.8419717
##
## Tuning parameter 'intercept' was held constant at a value of TRUE
##
# R2 es 0.84. Posee un buen ajuste.
# Graficamos los outliers
library(car)
plot(income, foodexp, pch=20, main = "Regresión Paramétrica")
abline(coef(lm(foodexp~income)), col='green')
out <- influence.measures(fitLS)
showLabels(engel$foodexp, engel$income, row.names(engel),
id.method=row.names(engel)[row.names(engel) %in% names(which(apply(out$is.inf, 1, any)))])
## 49 59 61 85 105 119 125 128 138 155 158 220
## 49 59 61 85 105 119 125 128 138 155 158 220
# Existencia de outliers. Lo conveniente es aplicar regresión robusta.
fitLSR <- ltsReg(income, foodexp)
summary(fitLSR)
##
## Call:
## ltsReg.default(x = income, y = foodexp)
##
## Residuals (from reweighted LS):
## Min 1Q Median 3Q Max
## -188.51 -36.66 0.00 34.26 161.44
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Intercept -5.25412 12.19031 -0.431 0.667
## income 0.68244 0.01261 54.133 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 69.49 on 209 degrees of freedom
## Multiple R-Squared: 0.9334, Adjusted R-squared: 0.9331
## F-statistic: 2930 on 1 and 209 DF, p-value: < 2.2e-16
Como el income(beta1), poseen un p-valor < 0.05 rechazamos la hipótesis nula por lo que solo él debe aparecer en la ecuación. El intercepto no pues aceptamos la hipótesis nula.
La ecuación quedaría como: foodexp = 0.68244 * income
F-statistic (ANOVA), el p-value es < 0.001 nos aporta que el modelo es significativo.
Se nos muestra en F-statistic del paso anterior. Seria: F(1,209)=2930, p-value < 0.001
La bondad de ajuste viene dada por R2 que en este caso es de 0.9334, el modelo es muy eficaz (entre 0 y 1).
# graficamos...
par(mfrow=c(1,2))
plot(income, foodexp, pch=20, main = "Regresión Paramétrica")
abline(coef(lm(foodexp~income)), col='green')
plot(foodexp~income, pch=20, main = "Regresión Robusta")
abline(fitLSR, col='red')
# Obtenemos los valores...
summary(fitLS)$r.squared
## [1] 0.8303646
summary(fitLS)$coefficients[,4]
## (Intercept) income
## 1.573614e-17 9.918947e-92
summary(fitLSR)$r.squared
## [1] 0.9334254
summary(fitLSR)$coefficients[,4]
## Intercept income
## 6.669065e-01 6.193203e-125
La regresión robusta obtuvo un mejor valor RSaquare, pasando del 0.83 a 0.93 con lo que ofrece un mayor ajuste.
Utilizaremos los datos “chicago”, del paquete “faraway”. Corresponden a datos de un estudio sobre la disponibilidad de seguros en Chicago desde diciembre 1977 a febrero 1978. Los datos están dados para cada código postal de Chicago (nombres de las filas o casos). Las variables son:
Queremos explicar la variable “involact” en función de las demás variables, excepto “volact”, y con “income” en logaritmo.
Nuevamente, en lugar de ajustar todos los modelos robustos que hemos visto, recomiendo utilizar la función “ltsReg” del paquete “rrcov” que permite obtener un resúmen del modelo comparable al de la función “lm”.
Comparar entre los modelos lm y ltsReg la significación de los coeficientes, el ajuste global y la bondad de ajuste.
# Obtenemos los valores...
library(faraway)
data("chicago")
head(chicago)
## race fire theft age volact involact income
## 60626 10.0 6.2 29 60.4 5.3 0.0 11744
## 60640 22.2 9.5 44 76.5 3.1 0.1 9323
## 60613 19.6 10.5 36 73.5 4.8 1.2 9948
## 60657 17.3 7.7 37 66.9 5.7 0.5 10656
## 60614 24.5 8.6 53 81.4 5.9 0.7 9730
## 60610 54.0 34.1 68 52.6 4.0 0.3 8231
attach(chicago)
# preparación de los datos. Aplicamos logaritmos sobre la columna income, eliminamos la columna no usada en la regresión.
datos <- chicago
datos$income <- log(datos$income)
datos <- datos[,-5]
head(datos)
## race fire theft age involact income
## 60626 10.0 6.2 29 60.4 0.0 9.371098
## 60640 22.2 9.5 44 76.5 0.1 9.140240
## 60613 19.6 10.5 36 73.5 1.2 9.205127
## 60657 17.3 7.7 37 66.9 0.5 9.273878
## 60614 24.5 8.6 53 81.4 0.7 9.182969
## 60610 54.0 34.1 68 52.6 0.3 9.015663
# graficamos a nivel general...
pairs(datos, pch=20)
# graficamos con correlación paramétrica
library(corrplot)
par(mfrow = c(1,1))
corrplot(cor(datos), type = "lower")
corrplot(cor(datos), method = "number", type = "lower")
# Gráfico resumen con distribucion (diagonal), linea de ajuste y scatter plot bivariable (inferior) y valor de correlacion y significancia (superior).
library(PerformanceAnalytics)
chart.Correlation(datos, histogram=TRUE, pch=19)
En el gráfico se aprecia una correlación fuerte y positiva entre involact y fire (0.70), involact y race (0.71) y poco significativa pero positiva entre involact y age (0.48)
Una correlacion fuerte y negativa entre involact e income (-0.70).
Como se aprecian outliers, realizamos la correlación robusta (en esta caso winsorizada).
library(WRS)
winall(datos)
## $cor
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1.0000000 0.7521896 0.5300843 0.2588940 0.7305137 -0.7818053
## [2,] 0.7521896 1.0000000 0.4285737 0.4674647 0.8580696 -0.8155316
## [3,] 0.5300843 0.4285737 1.0000000 0.3386472 0.4067537 -0.4810901
## [4,] 0.2588940 0.4674647 0.3386472 1.0000000 0.5294211 -0.5667558
## [5,] 0.7305137 0.8580696 0.4067537 0.5294211 1.0000000 -0.7360743
## [6,] -0.7818053 -0.8155316 -0.4810901 -0.5667558 -0.7360743 1.0000000
##
## $cov
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 678.375060 101.9276781 118.657817 91.087831 8.51273358 -3.25362950
## [2,] 101.927678 27.0682239 19.163367 32.853534 1.99736818 -0.67796175
## [3,] 118.657817 19.1633673 73.864015 39.315819 1.56406105 -0.66065901
## [4,] 91.087831 32.8535338 39.315819 182.476438 3.19970860 -1.22330262
## [5,] 8.512734 1.9973682 1.564061 3.199709 0.20017576 -0.05262134
## [6,] -3.253630 -0.6779618 -0.660659 -1.223303 -0.05262134 0.02553107
##
## $p.values
## [,1] [,2] [,3] [,4] [,5]
## [1,] NA 3.092819e-08 0.0002647333 0.0833663910 1.021741e-07
## [2,] 3.092819e-08 NA 0.0036608185 0.0014461527 1.162315e-11
## [3,] 2.647333e-04 3.660818e-03 NA 0.0228120339 5.934029e-03
## [4,] 8.336639e-02 1.446153e-03 0.0228120339 NA 2.699239e-04
## [5,] 1.021741e-07 1.162315e-11 0.0059340287 0.0002699239 NA
## [6,] 5.057372e-09 4.676761e-10 0.0010219429 0.0000857958 7.590220e-08
## [,6]
## [1,] 5.057372e-09
## [2,] 4.676761e-10
## [3,] 1.021943e-03
## [4,] 8.579580e-05
## [5,] 7.590220e-08
## [6,] NA
Se aprecian datos ligeramente superiores aunque siguen el mismo patrón. También en el caso de involact con income es fuerte y neamtiva (-0.73)
Se aprecia también mirando p-value, que en todos los casos, p-value < 0.05, por lo que rechazamos la hipótesis nula (correlación es 0). Existe una relación significativa entre las variables.
summary(fitLM0 <- lm(involact ~ ., data = datos))
##
## Call:
## lm(formula = involact ~ ., data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.85393 -0.16922 -0.03088 0.17890 0.81228
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.573976 3.857292 -0.927 0.359582
## race 0.009502 0.002490 3.817 0.000449 ***
## fire 0.039856 0.008766 4.547 4.76e-05 ***
## theft -0.010295 0.002818 -3.653 0.000728 ***
## age 0.008336 0.002744 3.038 0.004134 **
## income 0.345762 0.400123 0.864 0.392540
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3345 on 41 degrees of freedom
## Multiple R-squared: 0.7517, Adjusted R-squared: 0.7214
## F-statistic: 24.83 on 5 and 41 DF, p-value: 2.009e-11
El intercepto posee un p-value > 0.05 por lo que lo excluimos al no ser significativo. Volvemos a ajustar…
summary(fitLM <- lm(involact ~ .-1, data = datos))
##
## Call:
## lm(formula = involact ~ . - 1, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.87307 -0.14861 -0.01933 0.19550 0.81707
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## race 0.007994 0.001881 4.250 0.000116 ***
## fire 0.036445 0.007942 4.589 4e-05 ***
## theft -0.009563 0.002700 -3.541 0.000990 ***
## age 0.007065 0.002373 2.977 0.004810 **
## income -0.024709 0.015067 -1.640 0.108480
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.334 on 42 degrees of freedom
## Multiple R-squared: 0.8708, Adjusted R-squared: 0.8554
## F-statistic: 56.6 on 5 and 42 DF, p-value: < 2.2e-16
Ahora todas las variables son significativas, a excepción de income con p-value > 0.05. Pero le añadimos a la ecuación.
La ecuación quedaría como: involact = 0.0095 * race + 0.0398 * fire - 0.010 * theft + 0.008 * age - 1.64 * income
F-statistic (ANOVA), el p-value es < 0.001 nos aporta que el modelo es significativo.
Se nos muestra en F-statistic del paso anterior. Seria: F(2,42) = 101.2, p-value < 0.001
La bondad de ajuste viene dada por R2 que en este caso es de 0.87, el modelo es muy eficaz (entre 0 y 1). El modelo es significativo y bueno.
anova(fitLM0, fitLM)
## Analysis of Variance Table
##
## Model 1: involact ~ race + fire + theft + age + income
## Model 2: involact ~ (race + fire + theft + age + income) - 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 41 4.5882
## 2 42 4.6843 -1 -0.096073 0.8585 0.3596
# p-valor > 0.05 nos dice que eliminar el intercepto no tiene efectos significativos en el modelo.
avPlots(fitLM)
vif(fitLM)
## race fire theft age income
## 3.374141 6.258661 4.712647 9.820465 8.185995
# Como los datos que arroja la prueba no son mayores de 10, nos dice que los predictores NO están correlacionados. No hay colinealidad.
summary(step(fitLM, action = F, direction = "forward"))
## Start: AIC=-98.38
## involact ~ (race + fire + theft + age + income) - 1
##
## Call:
## lm(formula = involact ~ (race + fire + theft + age + income) -
## 1, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.87307 -0.14861 -0.01933 0.19550 0.81707
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## race 0.007994 0.001881 4.250 0.000116 ***
## fire 0.036445 0.007942 4.589 4e-05 ***
## theft -0.009563 0.002700 -3.541 0.000990 ***
## age 0.007065 0.002373 2.977 0.004810 **
## income -0.024709 0.015067 -1.640 0.108480
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.334 on 42 degrees of freedom
## Multiple R-squared: 0.8708, Adjusted R-squared: 0.8554
## F-statistic: 56.6 on 5 and 42 DF, p-value: < 2.2e-16
# En esta caso las variables se incluyen en el modelo final.
library(relaimpo)
(crlm <- calc.relimp(fitLM0, type = c("lmg", "last", "first", "pratt"), rela=T))
## Response variable: involact
## Total response variance: 0.4017299
## Analysis based on 47 observations
##
## 5 Regressors:
## race fire theft age income
## Proportion of variance explained by model: 75.17%
## Metrics are normalized to sum to 100% (rela=TRUE).
##
## Relative importance metrics:
##
## lmg last first pratt
## race 0.29243355 0.24877115 0.29140971 0.46388154
## fire 0.30950866 0.35299237 0.28272657 0.54707074
## theft 0.05749694 0.22790620 0.01280705 -0.07206704
## age 0.12207995 0.15757884 0.12945724 0.18789031
## income 0.21848090 0.01275145 0.28359943 -0.12677556
##
## Average coefficients for different model sizes:
##
## 1X 2Xs 3Xs 4Xs 5Xs
## race 0.013882353 0.010877318 0.009202480 0.008628879 0.009502223
## fire 0.047902503 0.040131301 0.035927286 0.035527385 0.039856040
## theft 0.004254602 -0.002893112 -0.005717252 -0.007625214 -0.010294505
## age 0.013356717 0.008219986 0.006461752 0.006660529 0.008335600
## income -1.798818702 -1.363378514 -0.889943051 -0.344752568 0.345761521
plot(crlm)
# En la mayoría de los métodos se aprecia como "race" y "fire" los que más importancia tienen y de ellos "fire" el más imporante. (que además son las de mayor significancia con ***)
par(mfrow = c(2,2))
plot(fitLM)
par(mfrow=c(1,1))
plot(fitLM, which=1)
# Están bastante cercanos a 0, no hay casi dispersion, aparecen 3 outliers.
plot(fitLM, which=2)
# Se aprecia que se ajustan lso residuos a la curva por lo que existe normalidad. Se muestran los 3 outliers
plot(fitLM, which=3)
# No se aprecian enfecto estilo embudo.
plot(fitLM, which=4)
# Se muestra la existencia de los 3 outliers.
outlierTest(fitLM)
##
## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferonni p
## 60610 -3.262772 0.0022278 0.1047
# Solo está detectando uno, el más alejado el 60610.
avPlots(fitLM)
influencePlot(fitLM, id.method = "noteworthy", main="Grafico Influencia")
## StudRes Hat CookD
## 60610 -3.2627716 0.2105705 0.46185190
## 60607 0.4509906 0.6383261 0.07318241
# Los más influyentes son los outliers 60610 y 60607 aunque con el outlierTest solo detectó el 60610.
qqPlot(fitLM)
# Ya comentado anteriormente que se ajusta bastante a la recta.
library(MASS)
sresid <- studres(fitLM)
hist(sresid, freq = F, main = "Distribucion de los residuos")
xfitLM <- seq(min(sresid), max(sresid), length=40)
yfitLM <- dnorm(xfitLM)
lines(xfitLM, yfitLM)
# Distribución bastante normal, aunque ligeramente sesgada a izda.
ncvTest(fitLM)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 9.205779 Df = 1 p = 0.002412523
# Se rechaza Ho por lo que no hay varianza constante en los errores.
spreadLevelPlot(fitLM)
##
## Suggested power transformation: 0.4330704
crPlots(fitLM)
# Race, Fire y age son lo que más se mantienen en la linealidad.
durbinWatsonTest(fitLM)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.1106017 2.210188 0.552
## Alternative hypothesis: rho != 0
# p-valor al límite 0.05.
pairs(datos)
out <- influence.measures(fitLM)
plot(fitLM, which=4)
# Existencia de outliers. Lo conveniente es aplicar regresión robusta.
# library(rrcov)
fitLMR0 <- ltsReg(involact ~ race + fire + theft + age + income)
summary(fitLMR0)
##
## Call:
## ltsReg.formula(formula = involact ~ race + fire + theft + age +
## income)
##
## Residuals (from reweighted LS):
## Min 1Q Median 3Q Max
## -0.387900 -0.077283 -0.005965 0.115227 0.496925
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Intercept -4.324e-01 3.313e-01 -1.305 0.2002
## race 5.678e-03 1.694e-03 3.351 0.0019 **
## fire 4.786e-02 5.861e-03 8.166 1.03e-09 ***
## theft -8.572e-03 1.919e-03 -4.467 7.57e-05 ***
## age 4.202e-03 1.914e-03 2.195 0.0347 *
## income 2.025e-05 2.115e-05 0.958 0.3445
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2202 on 36 degrees of freedom
## Multiple R-Squared: 0.8475, Adjusted R-squared: 0.8263
## F-statistic: 40 on 5 and 36 DF, p-value: 1.005e-13
# Todos menos income, poseen un p-valor < 0.05 rechazamos la hipótesis nula por lo que debeb aparecer en la ecuación. El intercepto no pues aceptamos la hipótesis nula.
fitLMR <- ltsReg(involact ~ race + fire + theft + age + income - 1)
summary(fitLMR)
##
## Call:
## ltsReg.formula(formula = involact ~ race + fire + theft + age +
## income - 1)
##
## Residuals (from reweighted LS):
## Min 1Q Median 3Q Max
## -0.39796 -0.04698 0.00000 0.10471 0.52330
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## race 4.580e-03 1.417e-03 3.233 0.002577 **
## fire 4.530e-02 5.858e-03 7.732 3.07e-09 ***
## theft -7.510e-03 1.899e-03 -3.955 0.000333 ***
## age 1.983e-03 1.463e-03 1.356 0.183368
## income -5.939e-06 6.326e-06 -0.939 0.353900
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2263 on 37 degrees of freedom
## Multiple R-Squared: 0.9072, Adjusted R-squared: 0.8947
## F-statistic: 72.38 on 5 and 37 DF, p-value: < 2.2e-16
Ahora todas las variables son significativas, a excepción de income con p-value > 0.05. Pero le añadimos a la ecuación.
Respecto al modelo paramétrico comentar que da menor significancia a la variable “race” e income no aparece con su estimado en netativo.
La ecuación quedaría como: involact = 0.00458 * race + 0.0453 * fire - 0.00751 * theft + 0.0019 * age - 0.0000059 * income
F-statistic (ANOVA), el p-value es < 0.001 nos aporta que el modelo es significativo.
Se nos muestra en F-statistic del paso anterior. Seria: F(5,37)=72.38, p-value < 0.001
La bondad de ajuste viene dada por R2 que en este caso es de 0.90, el modelo es muy eficaz (entre 0 y 1). El modelo es significativo y bueno.
Comentar que R2 el ligeramente inferior al obtenido por el anterior (paramétrico)
detach(chicago)
rm(list = ls())